#include <fstream>
#include <iostream>
#include <iomanip>
#include "TFile.h"
#include "TTree.h"
#include "TLine.h"
#include "TMath.h"
#include "TSystem.h"
#include "TStyle.h"
#include "TH1.h"
#include "TH2.h"
#include "TCanvas.h"
#include "TLegend.h"
#include "TEllipse.h"
#include "TClonesArray.h"


using namespace std;

void RPhi_75L(){

  int ret = gSystem->Load( "../Objects/libObjectClass.so" );
  if( ret!=0 ){
    cout << "Failed to Load Library: ret = " << ret << endl;
    return;
  }

  int zpos = -465;

  TH2F *hXZ[5];
  TH2F *hXY[5];
  TH1F *hPhi[15][5];
  for(int i=0;i<5;i++){
    hXZ[i]=new TH2F(Form("hXZ_%d",i),Form("hXZ_%d",i),300,-150,150,350,-650,-300);
    hXY[i]=new TH2F(Form("hXY_%d",i),Form("hXY_%d",i),150/2,-150,150,75/2,0,150);
    for(int j=0;j<15;j++)
      hPhi[j][i]=new TH1F(Form("hPhi_%d_%d",j,i),Form("hPhi_%d_%d",j,i),190,-5,185);
  }


  TFile *fread=new TFile( "/scratch/shinichi/OutputFiles/75L_RPhi/75L_Full_2phase.root" );
  TTree *tree=(TTree*)fread->Get("tree");

  int nPhoton,nPMThits[5];
  tree->SetBranchAddress("nPhoton",&nPhoton);
  tree->SetBranchAddress("nJPMThits",&nPMThits[0]);
  tree->SetBranchAddress("nMPMT1hits",&nPMThits[1]);
  tree->SetBranchAddress("nMPMT2hits",&nPMThits[2]);
  tree->SetBranchAddress("nMPMT3hits",&nPMThits[3]);
  tree->SetBranchAddress("nMPMT4hits",&nPMThits[4]);

  int nev = tree->GetEntries();
  cout<<"NEV = "<<nev<<endl;

  int nphi[15] = {3,5,8,11,14,17,20,24,27,30,33,36,39,19,19}; //305

  int NPHOTON=0;
  int NPMTHITS[5]={0};
  int IEV = 0;

  for(int iev=0;iev<nev;iev++){

    tree->GetEntry(iev);

    NPHOTON += nPhoton;
    for(int ipmt=0;ipmt<5;ipmt++) NPMTHITS[ipmt] += nPMThits[ipmt];
    
    if((iev%10)==9){

      int runid = IEV;
      int id=-1;
      int runID= runid % 305;
      if( runID>=0  && runID<nphi[0] )  id=0;
      else if( runID<(nphi[0]+nphi[1]) )     id=1;
      else if( runID<(nphi[0]+nphi[1]+nphi[2]) )  id=2;
      else if( runID<(nphi[0]+nphi[1]+nphi[2]+nphi[3]) )  id=3;
      else if( runID<(nphi[0]+nphi[1]+nphi[2]+nphi[3]+nphi[4]) )  id=4;
      else if( runID<(nphi[0]+nphi[1]+nphi[2]+nphi[3]+nphi[4]+nphi[5]) )  id=5;
      else if( runID<(nphi[0]+nphi[1]+nphi[2]+nphi[3]+nphi[4]+nphi[5]+nphi[6]) )  id=6;
      else if( runID<(nphi[0]+nphi[1]+nphi[2]+nphi[3]+nphi[4]+nphi[5]+nphi[6]+nphi[7]) )  id=7;
      else if( runID<(nphi[0]+nphi[1]+nphi[2]+nphi[3]+nphi[4]+nphi[5]+nphi[6]+nphi[7]+nphi[8]) )  id=8;
      else if( runID<(nphi[0]+nphi[1]+nphi[2]+nphi[3]+nphi[4]+nphi[5]+nphi[6]+nphi[7]+nphi[8]+nphi[9]) )  id=9;
      else if( runID<(nphi[0]+nphi[1]+nphi[2]+nphi[3]+nphi[4]+nphi[5]+nphi[6]+nphi[7]+nphi[8]+nphi[9]+nphi[10]) )  id=10;
      else if( runID<(nphi[0]+nphi[1]+nphi[2]+nphi[3]+nphi[4]+nphi[5]+nphi[6]+nphi[7]+nphi[8]+nphi[9]+nphi[10]+nphi[11]) )  id=11;
      else if( runID<(nphi[0]+nphi[1]+nphi[2]+nphi[3]+nphi[4]+nphi[5]+nphi[6]+nphi[7]+nphi[8]+nphi[9]+nphi[10]+nphi[11]+nphi[12]) )  id=12;
      else if( runID<(nphi[0]+nphi[1]+nphi[2]+nphi[3]+nphi[4]+nphi[5]+nphi[6]+nphi[7]+nphi[8]+nphi[9]+nphi[10]+nphi[11]+nphi[12]+nphi[13]) )  id=13;
      else if( runID<(nphi[0]+nphi[1]+nphi[2]+nphi[3]+nphi[4]+nphi[5]+nphi[6]+nphi[7]+nphi[8]+nphi[9]+nphi[10]+nphi[11]+nphi[12]+nphi[13]+nphi[14]) )  id=14;
      else if( id==-1 ){ cout<<runid<<endl;return;}
      
      int sum=0;
      if(id!=0) for(int ii=0;ii<id;ii++) sum+=nphi[ii];
      
      int iphi = (runID-sum) % nphi[id];
      int ir   = id;
      int iz   = runid / 305 ;
      
      //      cout<<NPMTHITS[0]<<" "<<NPHOTON<<endl;
      float rr = 5.+10.*ir;
      float phi = iphi*(TMath::Pi()/(float)(nphi[ir]-1)); //rad
      float zz;
      if(iz<24) zz = -595.+iz*10.;
      else      zz = -360.+iz*1.;

      if(zz==zpos){
	//	cout<<"("<<rr*cos(phi)<<","<<rr*sin(phi)<<")= "<<NPMTHITS[0]<<endl;
	for(int ipmt=0;ipmt<5;ipmt++){
	  hXY[ipmt]->Fill(rr*cos(phi),rr*sin(phi),NPMTHITS[ipmt]/(float)NPHOTON);
	  hPhi[ir][ipmt]->Fill(phi*180/TMath::Pi(),NPMTHITS[ipmt]/(float)NPHOTON);
	  hPhi[ir][ipmt]->SetBinError(hPhi[ir][ipmt]->FindBin(phi*180/TMath::Pi()),sqrt(NPMTHITS[ipmt])/(float)NPHOTON);
	}
      }
      if(iphi==0||iphi==nphi[ir]-1){
	for(int ipmt=0;ipmt<5;ipmt++)
	  hXZ[ipmt]->Fill(rr*cos(phi),zz,NPMTHITS[ipmt]/(float)NPHOTON);
      }
      NPHOTON = 0;
      for(int ipmt=0;ipmt<5;ipmt++) NPMTHITS[ipmt] = 0;
      IEV++;
    }

  }

  gStyle->SetOptStat(000000);
  gStyle->SetPadBottomMargin(0.1);
  //  gStyle->SetPadTopMargin(0.05);
  gStyle->SetPadRightMargin(0.12);
  gStyle->SetPadLeftMargin(0.1);

  int pos[5]={5,3,1,7,9};
  char cpmt[5][64]={"JPMT","Meg1","Meg2","Meg3","Meg4"};

  TCanvas *c1=new TCanvas("c1","c1",1200,1000);
  c1->Divide(3,3);
  for(int i=0;i<5;i++){
    c1->cd(pos[i]);
    char name[128];
    sprintf(name,"%s | Z = %d (mm)",cpmt[i],zpos);
    hXY[i]->SetMaximum(0.1);
    hXY[i]->SetTitle(name);
    hXY[i]->SetXTitle("X (mm)");
    hXY[i]->SetYTitle("Y (mm)");
    hXY[i]->Draw("colz");
  }

  TCanvas *c3=new TCanvas("c3","c3",1200,1000);
  c3->Divide(3,3);
  for(int i=0;i<5;i++){
    c3->cd(pos[i]);
    hXZ[i]->SetMaximum(0.1);
    hXZ[i]->SetXTitle("X (mm)");
    hXZ[i]->SetYTitle("Y (mm)");
    hXZ[i]->Draw("colz");
  }


  TCanvas *c2=new TCanvas("c2","c2",1200,1000);
  c2->Divide(3,3);

  TLegend *lg=new TLegend(0.2,0.1,0.9,0.9);
  lg->SetFillStyle(0);

  for(int i=0;i<5;i++){
    c2->cd(pos[i]);
    TH1 *frame=gPad->DrawFrame(-5,0,185,0.1);
    char name[128];
    sprintf(name,"%s | Z = %d (mm)",cpmt[i],zpos);
    frame->SetTitle(name);
    frame->SetXTitle("#phi (deg)");
    frame->SetYTitle("NEV");
    for(int j=0;j<15;j++){
      hPhi[j][i]->SetMarkerSize(1);
      hPhi[j][i]->SetMarkerStyle(21);
      hPhi[j][i]->SetMarkerColor(51+(int)(3.5*j));
      hPhi[j][i]->SetLineColor(51+(int)(3.5*j));
      if(j%2==1) continue;
      hPhi[j][i]->Draw("samepe");
      if(i==0) lg->AddEntry(hPhi[j][0],Form("R = %d (mm)",j*10+5),"lp");

    }
  }
  c2->cd(8);
  lg->Draw();

}

void write(){

  int ret = gSystem->Load( "../Objects/libObjectClass.so" );
  if( ret!=0 ){
    cout << "Failed to Load Library: ret = " << ret << endl;
    return;
  }


  TFile *fread=new TFile( "/scratch/shinichi/OutputFiles/75L_RPhi/75L_Full_2phase.root" );
  TTree *tree=(TTree*)fread->Get("tree");

  int nPhoton,nPMThits[5];
  tree->SetBranchAddress("nPhoton",&nPhoton);
  tree->SetBranchAddress("nJPMThits",&nPMThits[0]);
  tree->SetBranchAddress("nMPMT1hits",&nPMThits[1]);
  tree->SetBranchAddress("nMPMT2hits",&nPMThits[2]);
  tree->SetBranchAddress("nMPMT3hits",&nPMThits[3]);
  tree->SetBranchAddress("nMPMT4hits",&nPMThits[4]);

  int nev = tree->GetEntries();
  cout<<"NEV = "<<nev<<endl;

  float acc[5][305*34]={{0}};

  int    nphi[15] = {3,5,8,11,14,17,20,24, 27, 30, 33, 36, 39, 19, 19}; //305
  int sumnphi[15] = {0,3,8,16,27,41,58,78,102,129,159,192,228,267,286}; //305
  int NPHOTON=0;
  int NPMTHITS[5]={0};
  int IEV = 0;

  for(int iev=0;iev<nev;iev++){

    tree->GetEntry(iev);

    NPHOTON += nPhoton;
    for(int ipmt=0;ipmt<5;ipmt++) NPMTHITS[ipmt] += nPMThits[ipmt];
    
    if((iev%10)==9){

      int runid = IEV;
      int id=-1;
      int runID= runid % 305;
      if( runID>=0  && runID<sumnphi[1] )  id=0;
      else if( runID<sumnphi[2] )  id=1;
      else if( runID<sumnphi[3] )  id=2;
      else if( runID<sumnphi[4] )  id=3;
      else if( runID<sumnphi[5] )  id=4;
      else if( runID<sumnphi[6] )  id=5;
      else if( runID<sumnphi[7] )  id=6;
      else if( runID<sumnphi[8] )  id=7;
      else if( runID<sumnphi[9] )  id=8;
      else if( runID<sumnphi[10])  id=9;
      else if( runID<sumnphi[11])  id=10;
      else if( runID<sumnphi[12])  id=11;
      else if( runID<sumnphi[13])  id=12;
      else if( runID<sumnphi[14])  id=13;
      else if( runID<305 )         id=14;
      else if( id==-1 ){ cout<<runid<<endl;return;}
      
      int iphi = (runID-sumnphi[id]) % nphi[id];
      int ir   = id;
      int iz   = runid / 305 ;
      
      float rr = 5.+10.*ir;
      float phi = iphi*(TMath::Pi()/(float)(nphi[ir]-1)); //rad
      float zz;
      if( iz<24 ) zz = -595.+iz*10.;
      else        zz = -360.+iz*1.;

      for(int ipmt=0;ipmt<5;ipmt++)
	//acc[ipmt][IEV] = NPMTHITS[ipmt];	
	acc[ipmt][IEV] = NPMTHITS[ipmt]/(float)NPHOTON;

      NPHOTON = 0;
      for(int ipmt=0;ipmt<5;ipmt++) NPMTHITS[ipmt] = 0;
      IEV++;
    }
  }

  TH2F *test[5];
  TH2F *test1[5];
  TH2F *test2[5][10];
  for(int i=0;i<5;i++){
    test[i]=new TH2F(Form("htest_%d",i),Form("htest_%d",i),300,-150,150,300,-150,150);
    test1[i]=new TH2F(Form("htest1_%d",i),Form("htest1_%d",i),300,-150,150,320,-620,-300);
    for(int j=0;j<10;j++)
      test2[i][j]=new TH2F(Form("htest2_%d_%d",i,j),Form("htest2_%d_%d",i,j),300,-150,150,320,-620,-300);
  }
  cout<<"-----here-----"<<endl;
  
  float RR,theta;
  int ir,ir2;
  float xpos,ypos,zpos;

  for(int iz=0;iz<(-350-(-600));iz++){
    for(int iy=0;iy<150;iy++){
      for(int ix=0;ix<300;ix++){

	xpos = ix*1-150;
	ypos = iy*1;
	zpos = iz*1-600;
	RR   = sqrt(xpos*xpos+ypos*ypos);
	if( RR>=150.) continue;

	if( RR<=0. ) RR=0.01;
	theta = TMath::ACos( xpos/RR );

	//-----
	ir = (int)(RR/10.);
	if(ir>14){ cout<<ix<<" "<<iy<<" "<<iz<<endl;return;}

	float dphi = TMath::Pi()/(float)(nphi[ir]-1);
	float dphi2;
	if( ir==(int)(RR/10.-0.5) ){
	  if( ir==14 )  ir2 = ir;
	  else          ir2 = ir+1;
	  dphi2 = TMath::Pi()/(float)(nphi[ir2]-1);
	}
	else{
	  if( ir==0 )	ir2 = ir;
	  else          ir2 = ir-1;
	  dphi2 = TMath::Pi()/(float)(nphi[ir2]-1);
	}
	//-----
	int no1=-1;
	for(int j=0;j<nphi[ir];j++){
	  float tmpphi = j*dphi;
	  float delta = theta-tmpphi;
	  if( delta<0 ){ no1=j; break;}
	}
	if( no1<0 ) no1=nphi[ir]-1;

	int no2=-1;
	for(int j=0;j<nphi[ir2];j++){
	  float tmpphi = j*dphi2;
	  float delta = theta-tmpphi;
	  if( delta<0 ){ no2=j; break;}
	}
	if( no2<0 ) no2=nphi[ir2]-1;

	//-----
	int iza,izb;
	float dz,zza,zzb;
	if(zpos<-360){
	  iza = (int)((zpos+600)/10);
	  if(iza*10-595<zpos){
	    izb = iza+1;
	    if(iza==23) izb=iza;
	  }
	  else if(zpos<=-595)   izb = iza;
	  else                  izb = iza-1;  
	  dz = 10.*(iza-izb);
	  zza = -595.+iza*10.;
	  zzb = -595.+izb*10.;
	}
	else{
          iza = (int)(24+zpos+360);
	  if(iza==33)   izb = iza;
	  else          izb = iza+1;
	  dz = 1.;
	  zza = -360.+(iza-24)*1.;
	  zzb = -360.+(izb-24)*1.;
	}

	//-------
	for(int ipmt=0;ipmt<5;ipmt++){
	  
	  float Acc1_1a = acc[ipmt][305*iza+sumnphi[ir]+no1];
	  float Acc1_1b = acc[ipmt][305*izb+sumnphi[ir]+no1];
	  float Acc1_2a = acc[ipmt][305*iza+sumnphi[ir]+no1-1];
	  float Acc1_2b = acc[ipmt][305*izb+sumnphi[ir]+no1-1];
	  float Acc2_1a = acc[ipmt][305*iza+sumnphi[ir2]+no2];
	  float Acc2_1b = acc[ipmt][305*izb+sumnphi[ir2]+no2];
	  float Acc2_2a = acc[ipmt][305*iza+sumnphi[ir2]+no2-1];
	  float Acc2_2b = acc[ipmt][305*izb+sumnphi[ir2]+no2-1];

	  float Acc1a = (Acc1_1a-Acc1_2a)/ dphi*theta + (no1* dphi*Acc1_2a-(no1-1)* dphi*Acc1_1a)/dphi ;
	  float Acc1b = (Acc1_1b-Acc1_2b)/ dphi*theta + (no1* dphi*Acc1_2b-(no1-1)* dphi*Acc1_1b)/dphi ;
	  float Acc2a = (Acc2_1a-Acc2_2a)/dphi2*theta + (no2*dphi2*Acc2_2a-(no2-1)*dphi2*Acc2_1a)/dphi2 ;
	  float Acc2b = (Acc2_1b-Acc2_2b)/dphi2*theta + (no2*dphi2*Acc2_2b-(no2-1)*dphi2*Acc2_1b)/dphi2 ;
	  if(Acc1a<0) Acc1a=0;
	  if(Acc1b<0) Acc1b=0;
	  if(Acc2a<0) Acc2a=0;
	  if(Acc2b<0) Acc2b=0;

	  float Acca,Accb;	
	  if( ir==ir2 ){
	    Acca = Acc1a; 
	    Accb = Acc1b; 
	  }
	  else{
	    Acca = (Acc1a-Acc2a)/(10.*(ir-ir2))*RR + ((10.*ir+5.)*Acc2a-(10.*ir2+5.)*Acc1a)/(10.*(ir-ir2)) ;
	    Accb = (Acc1b-Acc2b)/(10.*(ir-ir2))*RR + ((10.*ir+5.)*Acc2b-(10.*ir2+5.)*Acc1b)/(10.*(ir-ir2)) ;
	  }
	  if(Acca<0) Acca=0;
	  if(Accb<0) Accb=0;

	  float Acc;
	  if(zza==zzb)    Acc = Acca;
	  else            Acc = (Acca-Accb)/(zza-zzb)*zpos + (zza*Accb-zzb*Acca)/(zza-zzb);
	  if(Acc<0) Acc=0;

	  //	  if(ipmt==0&&Acc<0) cout<<zza<<" "<<zzb<<" "<<Acca<<" "<<Accb<<" zpos="<<zpos<<" Acc="<<Acc<<endl;
	  if(ipmt==0&&ypos==0&&xpos==32&&zpos==-540) cout<<" Acc="<<Acc<<endl;

	  if(zpos==-351){
	    if( ipmt==0&&xpos==50.&&ypos==50. ){
	      cout<<"numz = "<<iza<<" "<<izb<<endl;
	      cout<<"R = "<<RR<<"   theta = "<<theta<<endl;
	      cout<<"1: "<<ir<<" " <<dphi*no1    <<" "<<Acc1_1a<<endl;
	      cout<<"2: "<<ir<<" " <<dphi*(no1-1)<<" "<<Acc1_2a<<endl;
	      cout<<"3: "<<ir2<<" "<<dphi2*no2   <<" "<<Acc2_1a<<endl;
	      cout<<"4: "<<ir2<<" "<<dphi2*(no2-1)<<" "<<Acc1_2a<<endl;
	      cout<<Acc1a<<" "<<Acc2a<<endl;
	      cout<<Acca<<endl;
	    }

	    test[ipmt]->Fill(xpos,ypos,Acc);
	    if(ypos!=0){
	      if(ipmt!=0)	  test[5-ipmt]->Fill(xpos,-ypos,Acc);
	      else          test[ipmt]  ->Fill(xpos,-ypos,Acc);
	    }

	  }

	  if(ypos==0) test1[ipmt]->Fill(xpos,zpos,Acc);

	  if(xpos>=3&&xpos<13){
	    test2[ipmt][(int)xpos-3]->Fill(ypos,zpos,Acc);
	    if(ypos!=0){
	      if(ipmt!=0)    test2[5-ipmt][(int)xpos-3]->Fill(-ypos,zpos,Acc);
	      else           test2[ipmt][(int)xpos-3]->Fill(-ypos,zpos,Acc);
	    }
	  }


	}
	
      }
    }
  }
	
		     
  gStyle->SetOptStat(000000);
  int pos[5]={5,3,1,7,9};
  char cpmt[5][64]={"JPMT","Meg1","Meg2","Meg3","Meg4"};

  TCanvas *c1=new TCanvas("c1","c1",1200,1000);
  c1->Divide(3,3);
  for(int i=0;i<5;i++){
    c1->cd(pos[i]);
    char name[128];
    sprintf(name,"%s | Z = %d (mm)",cpmt[i],-353);
    test[i]->SetMaximum(0.01);
    test[i]->SetTitle(name);
    test[i]->SetXTitle("X (mm)");
    test[i]->SetYTitle("Y (mm)");
    test[i]->Draw("colz");
  }

  TCanvas *c2=new TCanvas("c2","c2",1200,1000);
  c2->Divide(3,3);
  for(int i=0;i<5;i++){
    c2->cd(pos[i]);
    test2[i][7]->SetMaximum(0.01);
    test2[i][7]->SetXTitle("Y (mm)");
    test2[i][7]->SetYTitle("Z (mm)");
    test2[i][7]->Draw("colz");
  }

  TCanvas *c3=new TCanvas("c3","c3",1200,1000);
  c3->Divide(3,3);
  for(int i=0;i<5;i++){
    c3->cd(pos[i]);
    test1[i]->SetMaximum(0.01);
    test1[i]->SetXTitle("Y (mm)");
    test1[i]->SetYTitle("Z (mm)");
    test1[i]->Draw("colz");
  }


}

